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Abstract 


This paper is based on a previous paper [11] which introduced a methodology for extrapola- 
tion of time-series which contain time-varying frequency components. Here the methodology 
is more fully explained and additional examples of its use are presented. An important ad- 
dition is a comparison with published ARIMA extrapolation of the Air Passenger Data by 
SAS. See [16] for the SAS description of their results. 


The basis of the method is complex-demodulation [1] which decomposes a time sequence 
into its, potentially time varying, frequency components. These complex-valued time series 
show the variation of the original series in the neighborhood of a specific frequency. Thus, 
the method does not assume any form of stationarity of the process which generated the 
series. The method can be viewed as an exploratory method in order to reveal the nature of 
the variability in the time sequence. (See Tukey [18] for his original exposition of exploratory 
data analysis.) 


Inspection of the complex-demodulates will indicate whether the observed series could 
be treated as from a stationary process or whether there are significant changes in the 
behavior of the series over time and at some frequencies. For example, as is done here, the 
series may be demodulated at each of the “seasonal” frequencies in order to determine if 
there exists a predictable pattern in the seasonality. In addition, the relative amplitude of 
each demodulate indicates the importance of that frequency. 


The extrapolation method is as follows: 


1. The series is complex-demodulated at a selected set of frequencies. (See Bingham, 
Godfrey, and Tukey [1] for a description of complex-demodulation.) 


2. The amplitude and phase of each demodulate are extrapolated using linear extrapola- 
tion. 


3. The resulting extrapolated series are remodulated and summed to produce the extrap- 
olation of the original series. 


The use of linear extrapolation assumes that the time dependence of the series is such that 
the demodulates have smoothly varying amplitude and phase. However, each demodulate 
is treated separately. Therefore, each demodulate may have an assumed time-dependence 
which is independent of the others. 


A similar characterization of a non-stationary time-series was first given by Wald in [20] 
and [21] as part of his research concerning the seasonal variation of Austrian economic time- 
series. However, due to restricted computational resources, Wald assumed that there was a 
common time-dependence for the seasonal frequencies. The method is also closely related to 
Gabor’s time-frequency analysis, introduced in his 1946 paper, “Theory of Communication” 
[7] and further discussed in [8] and in the form of a lecture at MIT:[9]. 


Keywords: time-series, non-stationarity, complex-demodulation, seasonal variation, eco- 
nomic data, extrapolation. 
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Background 


The first version of this paper was presented at a Seminar at the University of Wisconsin, 
Madison, October 3-5, 1966 and published in Spectral Analysis of Time Series [11]. 


While complex-demodulation (See [1]) is becoming a widely-used method (sometimes 
also under the name “heterodyning” ) for the analysis of the time-frequency composition of 
time series, it appears not to have been pursued as a method for extrapolation. This is the 
main purpose of this paper. 


The original program implementation, discussed in [11], was written in Fortran for use 
on the Princeton 7094. It has been possible to re-compile the Fortran program used to 
carry out the computations described in the Example Applications Section. The data were 
also available except that the General Motors data which were only available in the original 
(graphical) form as I received them from the General Motors Research Center. These had 
to be “re-digitized.” 


Using these materials it was possible to recompute the reported results from 1966. 
These results are very close to what was originally reported, but differ slightly due to minor 
changes that I had made after the Conference, and likely due to changes in the system 
software since the 1960’s. After establishing this base, all of the Fortran was translated into 
Octave (Matlab) functions. These also matched the originally published results. 


However, it then became clear that a number of important points had not been fully 
thought through and implemented. The current implementation represents considerable 
improvement over the original. This code will be included as an Appendix and will be 
available on my web site when there is an opportunity to make the code better organized 
and more readable. 


Incidentally, the original programs were written in Fortran specifically for the Princeton 
University IBM 7094. This computing system had a capability which was unique at the 
time: a large CRT display and a 35mm film camera were attached online. And, most 
importantly this equipment was supported by excellent sortware written be Peter Warter. 
Thus, the machine could be used as a (rather large and somewhat cumbersome) graphical 
workstation, with high resolution graphical display and recording. This allowed interactive 
graphics-based computing that was not, at the time, available elsewhere. 


In the original computations the mean and trend were not removed!. Current results 
suggest that mean and trend removal improves the performance of the extrapolation. In 
principle, only the zero frequency demodulate would be affected by the trend, but, par- 
ticularly since the method is intended for relatively short series, leakage effects into the 
low frequency complex-demodulates can be substantial. The computations reported here 


! The trend and mean should be removed. I remember being aware of the leakage effects 
of trend, and remember providing for trend removal in related programs which I wrote 
at that time. Also, it may be appropriate to use trend removal and not include the zero- 
frequency demodulate in the extrapolation, or to consider the zero-frequency demodulate as 
a form of “trend.” It is also surprising that John Tukey did not remark on this point. I still 
have a draft with his comments (in red ball point, as usual), but they do not include mention 
of choices about trend removal. It is possible that he assumed I had “done the right thing.” 
However, the best choice about trend removal is not always clear. For the data reported 
here, trend removal appears to provide moderate improvement in the extrapolation. 
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include mean and trend removal. After the extrapolated series is computed, the mean and 
trend of the original is added back in. This procedure has made a substantial change to 
the low-frequency complex-demodulates and to the extrapolation. As in any analysis of 
systems which give rise to trends, the handling of trend-removal deserves careful thought 
and experimentation. When writing the 1967 paper I was aware of Wald’s work but un- 
aware of Gabor’s paper [7]. Now, Gabor’s work is better known and its connection with 
wavelet analysis has been established. I was also not aware of Ville’s fundamental paper 
(19]. (A translation in English is available, and so is the French original, reset using TRX 
and MetaPost.) 
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1. The Extrapolation Method 


1.1 Introduction 


This time series extrapolation method makes use of time-frequency localization. It is in- 
tended for use in situations where the change in the amplitude and phase of frequency 
components of a time series can be estimated and therefore can be extrapolated. 


One of the earliest developers of these ideas within the statistics community was Abra- 
ham Wald|[20][21]. Wald was then at the Austrian Institute for Business Cycle Research in 
Vienna, and was working on methods of seasonal adjustment of economic time series. His 
method, described in English in Godfrey and Karremann{12], is based on the assumption 
that the amplitude and phase of the fundamental seasonal frequency and its harmonics vary 
with time. He assumed that the variation within the set of seasonal frequencies derives from 
a single source. Use can be made of this correlation of the frequency components to reduce 
the complexity of estimation of the time varying seasonal pattern. This technique worked 
quite well in practice for economic data and was used by the Austrian Institute for some 
time from the mid 1930’s. The correlation of the frequencies was important to the practical 
feasibility of the method at a time before digital computers. This work was far advanced 
for its time and is still technically advanced when compared to methods in use today. In 
the late 1960’s the version presented in Godfrey and Karremann[12], was considered as a 
possible official method of seasonal adjustment by the US Bureau of Labor Statistics. Af- 
ter considerable review and experimentation, it was not adopted. One important negative 
factor was that the method tended to preserve more variation in the data than the widely 
used “Census Method.” Generally, short run variation is viewed by both producers and 
consumers of macroeconomic data, as a bad thing. However, this tends to lead to the un- 
fortunate fact that much of the information in the series is lost, particularly for the values 
near the end of the series. This is indicated by the fact that the magnitude of the monthly 
revisions of the data often exceeds the magnitude of monthly changes. Of course, there are 
other reasons for the large revisions, but the standard methods of seasonal adjustment tend 
to make this problem worse. 


A further, and better known, contribution was made by Dennis Gabor (1946)[7] in his 
1946 “Communication Theory” paper. He discussed time-frequency localization in the con- 
text of the analysis of auditory signals. Despite mentioning that the auditory mechanism 
naturally operated on a logarithmic scale, he used a linear scale for his frequency decom- 
position. Had he used a logarithmic decomposition the methods that he developed would 
have been, for all practical purposes, equivalent to a wavelet transform decomposition. 


Wald’s and Gabor’s ideas were developed before digital computer-based data analy- 
sis became practical. Gabor discussed and constructed analog computing systems for his 
analysis, but these were difficult to design, construct and operate. Once general-purpose 
digital computing was readily available, techniques, such as complex demodulation, became 
practical general-purpose tools. 


Many extrapolation techniques may be considered, based on the time-frequency local- 
ization assumption. Extrapolation of the complex demodulates is the one presented here. 
Another natural idea is to revisit Wald’s assumption by computing the eigenvalues and 
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eigenvectors of the n x m matrix of complex-demodulates, where n is the length of the series 
and m is the number of computed demodulates. It may be that most of the variation in 
the m vectors is contained in the first few eigenvectors. If this is the case, only estimates of 
these vectors are needed for extrapolation. 


While prediction theory and practice for stationary processes have been extensively 
developed (Wiener (1949) [23], and Whittle (1963) [22], applications of these methods have 
quite frequently been unsuccessful. In many of these cases it would appear that non- 
stationarity has contributed to the lack of predictive success. 


In this paper we will outline a practical procedure for forming predictions which requires 
only “time-local smoothness” of the generating process. In addition, this procedure will 
permit convenient treatment of processes with approximately periodic variations, such as 
certain seasonal patterns. This method also permits useful exploratory data analysis of 
the time series. Each complex-demodulate shows how the amplitude and relative phase of 
the series changes, at each selected frequency, as a function of time. This can not only 
indicate “predictable” changes, but can locate time-points where the structure of series 
changes suddenly. Such changes may indicate that the entire series should not be treated as 
having smoothly varying time-frequency structure, but should be sub-divided into separate 
sections. 


1.2 Prediction Model 


The usual representation of a time-dependent random sequence is:! 


tt = / etk(w, t)dw, (124) 
where k(w,t) is a random function of w and t which is orthogonal on w, and 2; is the time 
series which will be predicted from a sample realization. Thus, the time-dependent spectrum 
which is to be estimated takes the form 


f(w,t) = E(k(w, t)k* (w, t)), (1.2.2) 


where * denotes the complex conjugate. 


If the process were stationary, such that f(w,t) = f(w), an estimate, f(w), of the spectral 
density 
f(w) = E(k(w)k*(w)) (1.2.3) 


would be derived using standard methods. Since f(w) is not a function of time, there 
is no problem in using it in order to make inferences about future values of the series. 
However, in the non-stationary case considered here, the spectral density f(w,t) changes 
within each realization. Therefore, it is not reasonable to estimate the spectral density by 
taking a simple average down the realization since this results in an average which treats 
each observation as if it came from a time-independent system. 


As a simple example, consider the case of a linear trend in the frequency components 
of the spectral density function. Then the spectral density is: 


f(w, t) = g(w) + A(w)t. (1.2.4) 


' An alternative representation has been given by Priestley[15]. 
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Averaging down a realization of n observations from t = 1,...,n, yields the density function 


Hie eat * a (1.2.5) 


This is not helpful for the purpose of extrapolation of the series. 


If we are interested in prediction of x, at some point in time n +7, then an estimate 
of the spectral density of the process up to the time n+ 7 is required. In the simple case 
of a linear trend in the spectrum, an estimate based on (1.2.5) will be closest to the true 
changing spectrum at t = “ and farthest at t = 1 and t = n. This is clearly an undesirable 
property for prediction. 


The traditional procedure for estimation of spectra, and subsequent computation of 
optimal predicting filters has been based on estimation of serial correlation coefficients. An 
assumption of stationarity over the entire length of the realization is required in order to 
make the usual inferences from these estimates. An alternative procedure for estimation of 
spectra is to apply digital filtering. In any case, both these methods require stationarity. 


The spectrum of a process may be estimated by first band-pass filtering the realiza- 
tion at a set of frequencies (often, but not necessarily, equispaced) on the range 0 to 7, 
and then computing the variance of each filtered series. The band-pass filtering does not 
require stationarity, but computing the variance over the entire sample does. For the non- 
stationary case we replace the simple variance with a time-dependent series computed for 
each frequency. 


1.3 The Prediction Algorithm 


In making use of digital filtering, it seems natural to consider extrapolating the results of 
the digital filtering at each frequency. While the best procedure will very likely depend to 
some extent on particular properties of the series being studied, a fairly general procedure 
is to use an expression of the form: 


A(w,t) = K(w) + a(w)t + S>as(w)A(w,t — 8), (1.3.1) 


where K(w), a(w), and a,(w) are constants, and A(w,t) is an estimate of the power spectrum 
at frequency w and time ¢. A similar expression is used for the phase of the filtered series. 
A computationally simple way of doing this is as follows: 


First complex-demodulate? the series x, by forming 
xi4= L(e™i*z,), P= 0,220,735 t= Leen; (1.3.2) 


where «,; is the complex-demodulate of x; at frequency w; and L(_) is a low-pass filter. The 
m-+1 frequencies which are demodulated may, if we have no prior knowledge about the series, 
simply be chosen to lie uniformly between 0 and 7. The appropriate value of m will depend 
on estimation criteria to be discussed later. However, if we have some knowledge of the 
structure of the generating process of the realization, the demodulation frequencies should 


? Presentations of complex-demodulation may be found in Tukey(1959)[17], Bingham, 
Godfrey, Tukey(1967)[1] and Godfrey(1965) [10]. 
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be chosen in order to obtain estimates centered on frequencies of particular importance. 
In particular, if there seems to exist a periodic variation with a changing pattern, then 
demodulation should be done at the fundamental frequency and each of its harmonics. If 
the intent is to produce a result which does not include these quasi-periodic effects (as is the 
case in seasonal adjustment) the complex-demodulates at these frequencies can be excluded 
when computing the extrapolation. 


Based on (1.2.2) observe that 


f(wj,t) = E(2j,225 ,). (1.3.3) 


jt 


Thus, estimates of the changing spectrum may be formed by taking time-local moving av- 
erages of x;,a;,. Then, extrapolated values of x, are formed by extrapolating each complex- 
demodulate and remodulating and summing the series. The z,;; are extrapolated based on 
(1.3.1), by computing: 
a5¢ = |X; ,y| 
and 
jt = tan (—S[2;,2], R[x;,2)) 


using 


ly Ig 
Qj,t = K(wj) + a(wy)t t+ > as(ws)bjt-2 + > Bs (wz) bj,2— 
s=1 


s=0 
ly Ig 

|@5,el? = K (wy) + a(wy)t + D> ors(wy)|Bj,e-01? + >> Bo(wy)laj,e-01°- (1.3.4) 
s=1 s=0 


Carry out the extrapolation of the amplitude and phase to t = n+7, where 7 is the number 
of time periods into the future which are to be forecast and n is the time point of the last 
observation. Then remodulate the @; by 


je = Fie (1.3.5) 


and form the final prediction by 
a, = (3-8), (1.3.6) 
j=0 


where () is the real part of its argument. Since by (1.3.4) the variance is estimated within 
each local time interval, nothing in the procedure requires stationarity of either first or 
second moments. All that is required for success in prediction is predictable change in the 
complex-demodulates. If the original series is from a stationary white noise process then 
this procedure not prove more successful than the simple procedure of predicting the mean 
of the series. However, the more time varying structure there is in the process, the better 
will be the forecasts by this procedure. In particular, if there is a deterministic component 
given by 

x, = (a+ Bt)e 5", (137) 


then its complex-demodulate at w,; will be 
vi, = at pt. (1.3.8) 


Given a long enough realization to permit the formation of reasonable estimates of a and 
8, prediction of this component can be done with a very small error. 
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There are surely cases such that some frequencies are more important than others It is 
easy to introduce a weighting function which represents the relative importance of errors at 
each frequency. To implement this, equation (1.3.6) is replaced by 


a, = (Sees), (1.3.9) 
j=0 


where the c; are weights constrained by 


xe =1. (1.3.10) 


j=0 


The m+ 1 weights should be chosen to best reflect the relative costs of prediction errors at 
each frequency. 


1.4 Details of the Estimation Procedure 


In order to form predictions using the technique outlined above, choices must be made 
concerning the bandwidth of the low-pass filter used in demodulation and the coefficients 
for the extrapolation of the demodulates. These choices are clearly interdependent. For 
example, the bandwidth of the filter determines the number of observations which are lost 
at the end of the demodulates and their apparent smoothness. In any case, the choices 
depend on the behavior of the time series being considered and on the intended uses of the 
predictions. 


The basic procedure indicated above did not specify how to handle the end terms of 
the series when applying the low-pass filter and the extrapolation shown in (1.3.4). 


There are various ways to handle the original and partially filtered observations near 
the end of the realization. It appears that the most important consideration is to avoid 
introducing artificial changes to the phase of the resulting complex demodulates. The 
specific choice, described below, was arrived at after considerable experimentation. In the 
end what appears to be a quite robust set of parameter values was found. These seem to 
work over a wide range of signals. But, of course, they will not be appropriate for all types 
of non-stationarity. 


The computational steps for the algorithm are as follows: 


1. The complex-demodulates are formed from convolutions of unweighted moving averages. 
This generates the required low-pass filter using many fewer numerical operations than 
would be required if the weights given by the Fourier transform of the filter function were 
used. In addition, this partially (after some tapering) preserves the end terms at each 
step of the convolution for use in the prediction equation (1.4.7). Thus the following 
computations are implemented: 


to,52 = ze", $= 1, 2ie5 j=0,...,m 
k 
1 
C156 = Gren, De Bagitte t=k+1,...,n—k 
j=9,. »™m 
al eee 
(1.4.1) 
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If r convolutions of length 2k +1 are applied using (1.4.1), a series of n terms remains 
which has at each end rk terms which are not fully smoothed. The first and last & terms 
were not smoothed at all, the terms that are from 4 +1 to 2k and n—(k+1) to n— 2k 
were only convolved once, and so on to rk and n— rk. 


Then the first and last kr terms are replaced by 


Ur jt = Lr jrk+1 t=1,...,rk 
Lr jt = r,j,m—rk t=m—rk+1,...,m. 
(1.4.2) 


This has the effect of reducing the variability due to the end terms while also reducing 
modification of the phase. 


There are other possibly reasonable choices about the handling of the end terms of the 
series. For example, the summation in (1.4.1) could always be over the full length of 
the series. The most important consideration, it appears, is to ensure that the phase of 
the resulting series is not distorted in the end terms just before the extrapolation starts. 
This is an important area where further work is needed. 


2. Next the set of m+1 complex valued series x,,; for 7 = 0,1,...,m, must be extrapolated 
to a specified time point n+7. Since the series are necessarily fairly smooth functions of 
time it is natural to apply a low-order linear regressive filter. However, before applying 
the linear regressive filter, the trend coefficients are computed according to: 


ijt = Ay + bt + Yy.t- (1.4.3) 


The trend removed series is then used in the extrapolation. Finally, after the extrapo- 
lation, the trend is added back and the predicted complex-demodulates are computed 
by: 


ly lo 
X54 — S- Aslj,t—s + PS Bs41Yj,t—s> (1.4.4) 
s=1 s=0 
where y;, = 0 for t > n and then 
Bit = X54 + aj + byt. 


Since it is the magnitude and phase of x) which will be slowly changing functions of 
time, equation (1.4.4) should be interpreted in terms of polar coordinates. Thus it is 
appropriate to extrapolate |x;1|* and ¢;4 rather than the real and imaginary parts of 
x;- ?;,4 iS, of course, given by 


5, = tan—'(S[a;,2], R[x;,t]) (1.4.5) 


where tan~! is a function of two arguments which are respectively the imaginary and real 
ordinates of the angle ¢;;. In order to carry out the extrapolation, the computed phase 
must be “unwrapped” so that it does not switch from —z to 7. This must be done by 
choosing a tolerance value to decide if the signal has really crossed the +m boundary and 
therefore should have +27 added to it. Currently, this tolerance is set to 1.47. Smaller 
values of the tolerance value seem to cause a “ringing” effect. 


At this point it is necessary to choose values for the coefficients 11, 12, a, and 6,. These 
determine the properties of the transfer function given by (1.4.6). The a, and #, should 
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be chosen in order to pass the relatively low frequency variation in the amplitude and 
phase of the signal at each frequency. The transfer function is given by 


. Bsz°rA(n + 8 —t) 


T(z) = i, (1.4.6) 
1 = 7 G2? 
where 
Al) =1 1>0 
= 1<0. 


The behavior of T(z) should also be chosen depending on the number of time periods 
which are to be predicted. If very long-term forecasts are of most value then the transfer 
function should essentially be a low-pass filter. 


The a, and 8, could be estimated from the complex-demodulates using regression pro- 
cedures. However, in the present study we have simply chosen coefficients on the basis 
of general gain and phase criteria. 


3. Once the series @;; have been formed for values of t to n+ 7, where 7 is the number of 
periods which are to be predicted, then 


a, = Tey ages (1.4.7) 
j=l 


is, after applying the previously removed trend, the forecast series. However, there is one 
final adjustment which tends to reduce the overall bias in the forecast. This is based on 
the fact that, in principle, the first forecast value should be very close to the last observed 
value. Therefore the difference between these two values is computed as an offset and 
subtracted from the extrapolated series. Note that the summation is from j = 1 since 
the zero frequency complex-demodulate should not enter into the extrapolation as it is 
taken care of by trend removal. 


We have experimented with a number of parameter values. For the airline passenger 
data reported later r = 1 and kj = 11 were used. Other choices were tried, but they did 
not make much difference. It appears that, if there are strong harmonic components, which 
may be changing with time, the extrapolated values depend only weakly on the smoothing 
parameters. The method easily extrapolates the nearly periodic components in any case. 
Conversely, if there are no substantial nearly periodic components the extrapolation quickly 
diverges from the actual values regardless of the choice of smoothing parameters. 
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2. Example Applications 


This Section shows the results of the method for 2 datasets: (1) the frequently used “Inter- 
national Air Passengers” data, and (2) Buick Sales data used in the original paper[11]. The 
results for the Buick Sales Data, as originally reported in{11], are summarized in Section 4. 


2.1 Air Passenger Data 


When the materials for this paper were found, a time series of the numbers of international 
air passengers was also found. This was a monthly series of airline passenger numbers (in 
thousands) from January 1949 through December 1960. These data were first reported in 
Brown|/4], and have become a standard example in texts on time series analysis. The series 
actually dates from January 1946, so I recovered the first 3 years from the CAA Statistical 
Handbook of Civil Aviation[5]. 


The Air Passengers data were published in: 


1. CAA Statistical Handbook of Civil Aviation, Civil Aeronautics Administration, 1948-58. 
(Volumes for 1951-52 are missing from the Stanford Library.) The name until sometime 
between 1959 and 1961 was Statistical Handbook of Civil Aviation|5], 


2. FAA Statistical Handbook of Aviation|6], Federal Aviation Agency, 1961-62. This publi- 
cation continued after 1962, but no longer reported monthly air passenger data. Because 
of the lag in publication, the series ends with December 1960. 


As is usual with economic data, the series were revised in succeeding publications. 
Brown|4] first studied these data. He used the un-revised data, and his published dataset 
seems to have been widely copied in the literature. The revised data are used here. 


2.1.1 Data Analysis 


Based on the description above and details given in Appendix 1 (p. 34), the series from Jan. 
1946 through Dec. 1960 is composed of 180 monthly observations. The series has an obvious 
seasonal component. It is much more “regular” than the Buick Sales data, but still obviously 
not stationary. I do not recall trying these data with the original program. First, we show a 
set of three extrapolations, starting respectively at Jan. 1952, Jan. 1956, and Jan. 1959. All 
results were computed after trend removal. The frequencies used were m = 0,...,19. Since 
the data are sampled monthly, the set of seasonal frequencies is m = 1,...,6. Experimenting 
with m sequences with upper values of 12, 18, 24, and 30 showed that the value 24 (25 
frequencies including 0) gave the lowest RMS extrapolation error. This led to a spacing of 
the seasonal demodulates as every forth demodulate. This effect was not observed for the 
Buick Sales data, which was sampled at 36 samples per year. The higher sampling rate 
may explain why additional demodulates were not useful for the Buick data. Or the cause 
could be the irregularity of the series. The value of k for the filtering of the demodulates 
was 5 and the number of convolutions, r, was 1. 


The smoothing equation used was the same as for the Buick Sales data: 


£54 = 1.5725 4-1 _ 0.725 1-2 +0.07y%-1+0.06y-2, t=4k+4+2, ton, (2.1.1) 
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X54 = 1.5725 4-1 — 0.7%; 4-2, t=nt+ 1, to n+ T; (2.1.2) 


The following Table shows the series lengths and the extrapolations that were used. 


Case Length Extrapolation k r 
start end 

extr_1953 180 84 156 11 1 

extr_1956 180 120 120 11 1 

extr_1958 180 144 96 11 1 


Table 2.1.1: Table of Parameter Values. 


The end dates for the extrapolations were chosen so that they all end on the same date. 
This facilitates comparison of the plots. This series is easier to extrapolate than the Buick 
Sales data, except for the reduction in the long-term growth rate. The year 1958 was chosen 
to show one of the limits to extrapolation of short series. The same procedure was tried, but 
using a start date of 1949 (as was used in many other studies). This substantially reduced 
the upward bias in the extrapolation since rapid growth only started in about 1949. (TWA 
and PanAm started world-wide services at this time, using Lockheed Constellations and a 
bit later Douglas DC-6’s, then followed in 1958 by the Boeing 707.) 


There are two other interesting points about these three extrapolations: 


1. The extrapolation from 1953 fairly quickly (during 1956) loses the seasonal pattern 
from the first seven years. This could be due to the irregularity of the data during the 
first 2 years. The extrapolations from 1956 and 1958 do not show this behavior. 


2. The extrapolation from 1958 is somewhat worse than the one from 1956. This raises 
a general point: In non-stationary systems, more data may not just do little good, it 
can do harm. It appears that irrelevant data from the past had a negative influence in 
the 1958 extrapolation. This suggests the possible usefulness of assigning less weight 
to data from farther in the past, or just moving the start point ahead, so that a fixed 
time-window is used. 
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Figure 2.1.1: Original Series 
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Figure 2.1.2: Original Series log transform and Trend Removed 
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Figure 2.1.3: Complex-demodulates (4;,.) to t = 84 (1953) 
Series after trend Restored 


Note that in this Figure the phase plot for 7/2 was set to zero after 1960 due to the fact 
that the amplitude was reduced to zero at that point and thus would make no contribution 
to the forecast in any case. 
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Figure 2.1.4: Remodulated (¢ < 1953) and Predicted (t > 1953) Values 
Series after trend Restored 
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Figure 2.1.5: Difference: Original — Remodulated Series with Trend Removal 


DRAFT: 4.1 15 September 2022. 


Example Applications 13 


Next we consider the extrapolation of this series starting the extrapolation at 1956. 
The following 3 plots show the results. 
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Figure 2.1.6: Complex-demodulates (#;,,) to t = 120 (1956) 
Series after trend Restored 
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Figure 2.1.7: Remodulated (¢ < 1956) and Predicted (t > 1956) Values 
Series after trend Restored 
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Figure 2.1.8: Difference: Original — Remodulated Series 
Series with Trend Removal 


DRAFT: 4.1 15 September 2022. 


Example Applications 15 


Next we consider the extrapolation starting at 1958. The following 3 plots show the 
results. 
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Figure 2.1.9: Complex-demodulates (#;,) to t = 144 (1958) 
Series after trend Restored 
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Figure 2.1.10: Remodulated (¢ < 1958) and Predicted (t > 1958) Values 
Series after trend Restored 
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Figure 2.1.11: Difference: Original — Remodulated Series 
Series with Trend Removal 
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2.1.2 Comparison with Published Extrapolations 


The SAS website [16] shows the extrapolation of the monthly Air Passenger data to March 
1963. They started their computation at January 1949. They use the Box-Jenkins ARIMA 
program. (The Box-Jenkins method is described in [2] and on Wikipedia under Box- 
Jenkins). 


The company Spider Financial has also published their extrapolation using the Box- 
Jenkins ARIMA model on their web site. Their extrapolation of the Air Passenger data 
is from the end of the Series at January 1961 to March 1963. They state that they based 
their method on the Box-Jenkins method as described in Time Series: Forecast and Control 
by Box, Jenkins and Reinsel [3]. It appears that they applied the standard Box-Jenkins 
technique which also appears to be the same ARIMA model as used by SAS. 


First we show the remodulated and predicted values using complex-demodulation for 
the prediction starting at January 1961 and extending to March 1963 as in the Spider 
Financial extrapolation. Our data start at January 1946 whereas the Spider data start at 
January 1949. The difference in the start dates does not appear to be significant. 
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Figure 2.1.12: Complex-demodulates (4;) : t = 1946 to March 1963 


DRAFT: 4.1 15 September 2022. 


18 Prediction for Non-Stationary Stochastic Processes — II (2022) 


The following Figure shows the Spider Financial (SAS ARIMA) result along with the 
original, remodulated and extrapolated series. 
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Figure 2.1.13: Remodulated (t < 1961) and predicted (¢ > 1961) values (red) 
and Spider Financial predicted (green). Series after trend Restored 


It is clear, without showing the confidence bounds computed by Spider Financial or 
SAS, that there is practically no difference between the two ARIMA extrapolations and ours. 
However, in our case the extrapolation was done after choosing the obvious demodulation 
frequencies and just one parameter (the window length nl). No parameter fitting was 
required. Simple experiments with ni show that the extrapolation is only very slightly 
changed for values within the range 6 < nl < 14. Actually, nl = 8 is slightly closer to the 
Spider Financial result than nl = 11. The difference is certainly not significant in any sense. 
In addition, the demodulation method can be applied unchanged as new data are made 
available. 
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2.2 Buick Sales Data 


The “Buick Sales Data” were obtained from the General Motors Research Center, Warren, 
Michigan. They represent sales in units, at 10-day intervals, of Buick cars from 1961 through 
1965. Included late in this period was a major strike which halted production. 


This series was chosen partly due to its apparent time-varying seasonal pattern. It was 
also chosen due to its highly irregular behavior. 


2.2.1 Data Analysis 


The sampling interval was described as “approximately 10 day intervals.” However, data 
were recorded in such a way that there are exactly 36 observations per year. Thus the 


fundamental seasonal frequency and its harmonics occur at: 
kr 

= — cou LS. 2.2.1 

eT yee (2.2.1) 


We used these 18 frequencies and w = 0 for the demodulation frequencies. We present a 
total of three extrapolations, extrapolating from 1963, 1964, and 1965. The following Table 
shows the values of the parameters: 


Table 2.2.1: Table of Coefficient Values. 


Case Length Extrapolation k r 
start end 

extr_1963 159 72 108 2 1 

extr_1964 159 108 72 5 1 

extr_1965 159 144 36 2 1 


The values of k and r are shown above. The values of the extrapolation coefficients were: 
a = [1.57,—.7]; 8 = [.0,.07,.06]. (Therefore, in (1.4.4) 1, and lz were both 2.) These were 
arrived at by experiments and computation of the relevant transfer functions. Figure 2.2.3 
shows the demodulates and extrapolations of the demodulates for three selected frequencies. 
The graphs show the smoothed demodulates up to t = 144 and the extrapolated series after 
that point. The equation used was: 


X54 — 1.5725 4-1 = 0.725 4-2 + 0.07y:~-1 + 0.06yr~-2, (2.2.2) 


for t=4k+ 2, to n, and 
£5,t = 1.57%54-1 — 0.725,t-2, (2.2.3) 


fort=n-+1, ton+r. 


For the results shown here the values of n and 7 are: n = 144, r = 36. (7 = 20 was used in 
the original paper.) The length of the series is 159. Using our new notation the filter length 
was 2, and number of convolutions was 1. The filter length of 2 is as close as possible to the 
old filter length of 6. 2 in the new definition would be 5 in the old, total length, definition. 
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Figure 2.2.1 shows the original series, Figure 2.2.2 shows the original series from which 
the demodulates were computed. Figure 2.2.3 shows three of the demodulates of the series. 
Figure 2.2.4 shows the remodulation of the smoothed and extrapolated demodulates (#;,1). 
Finally, Figure 2.2.5 shows the difference between the original series and #;. We had 159 
observations of the original series z,. However, we used only the first 144 observations in 
computing #,;. Therefore, there are 15 points for which we have known values of 2; for 
comparison with the predicted values. The point at which we stopped using the original 
series for the computation is indicated by a vertical line in each figure. 


17-Nov-2020: gmdata_1961_66_5m_extr_1963 


10000 - 


8000 |- 
s 6000} 
=] 
a 
= | 
(= 
= 
© 4000 + 
2000 |. 
oO 1 1 1 1 1 
1961 1962 1963 1964 1965 1966 
Time 
Figure 2.2.1: Original Series 
17-Nov-2020: gmdata_1961_66_5m_extr_1963 
4400 - 
oa 
a2 
=> 
2 
= 2200 + 
c 
— 
os 
D 
= 
a=) 
i— 
Cc 
Zo Or 
iJ 
= 
= 
s 
xs 
= 
© -2200 + 
> 
fo) 
4400 1 1 1 1 \ 
1961 1962 1963 1964 1965 1966 


Time 


Figure 2.2.2: Original Series: Trend Removed 
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Figure 2.2.4: Remodulated (t < 72) and Predicted (¢ > 72 (1963) Values 
Series after Trend Restored 
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Figure 2.2.5: Difference: Original — Remodulated Series 
Series after Trend Restored 
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Figure 2.2.6: Complex-demodulates (#;,,) to t = 108 (1964) 
Series after trend Restored 
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Figure 2.2.7: Remodulated (t < 108) and Predicted (¢ > 108) Values 
Series after Trend Restored 
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Figure 2.2.8: Difference: Original — Remodulated Series 
Series after Trend Restored 
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Series after trend Restored 
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Figure 2.2.10: Remodulated (t < 144) and Predicted (t > 144) Values 
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Figure 2.2.11: Difference: Original — Remodulated Series 
Series after Trend Restored 


The accuracy of the predicted values from ¢ = 145 through ¢ = 159 (in Figure 2.2.4 
this is from the first vertical blue line to the end of the original series) is fairly similar to 
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what was shown in the first publication. In that paper the extrapolation was done for 20 
time points. For Figure 2.2.2 above, the extrapolation was computed for 36 points. This 
additional extrapolation shows that the seasonal variation is substantially retained in the 
extrapolated series. But, the seasonal peak values are under-estimated. This change in 
behavior appears to have occurred after about 1958. 


Note also that the plotted phases of the complex demodulates are different in the new 
computation. This is due to the fact that the old computation plotted the phase before the 
“unwinding” had been done. This was not intended. 


There is some evidence that the use of the series without trend removal permitted 
extrapolation of frequencies below the frequency with period 7 = 144. But, this may just be 
an artifact of this particular data set. 


2.3 Use of the Octave Scripts and Functions 


The code that computes these results is currently written in Octave. This is the most 
suitable and most widely available system for analysis of this kind. 


In very brief outline, here are the scripts that were used: 


1. nonstat_pred: This computes all of the results required to form the extrapolation. These 
results are stored in files so that they may be used by the various plotting scripts. When 
the script is run is requests a file. This file should contain, in Octave language form, 
all of the parameters and data that will be needed. An example file is included with 
the scripts and functions with name: inpt.m. Thus, 


octave:1>nonstat_pred 
octave:2>inpt.m 


will initiate processing of the data contained in inpt.m. The code prompts for any 
changes to the parameters. If no changes are required, “cr” will start the processing. 
All results are written to files so that display of the outputs may be carried out within 
the current Octave session, or may be continued later simply by entering Octave (in 
the same directory). 


2. A series of scripts whose names all start with plot_ may be used to show the results. 
Thus, plot_orig displays the original time series, plot_demods displays selected complex 
demodulates, and plot_orig norm displays the original series and the remodulated and 
extrapolated series. These are the scripts which were used for the plots shown here. 
The script all_plots generates PDFs of all the plots as shown in Section 2 (p. 8). 


For more details it may be best for now to look at the Octave source code. 


2.4 Further Developments 


It is important to develop more formal tests of the practical performance of this procedure. 
However, the classical tests of prediction error cannot be directly applied as they rely on 
stationarity. Since non-stationarity is just defined as the lack of stationarity it will take 
some work to find more specific definitions for which tests could be implemented. 


There are also a number of aspects of this method which deserve more study: 
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1. There is a trade-off between recent and more distant structure in the data. At each 
frequency, it takes a number of cycles to estimate the behavior at that frequency. This 
means that lower frequency variation must depend on more distant past. But, just 
using the frequency ratios to determine the distance is not likely a best choice. 


2. It is possible in some cases that the complex-demodulates are correlated. This was an 
assumption in Wald’s analysis. This can be tested, and used, by application of principal 
components methods. In an extreme case where one factor controls the behavior of all 
of the harmonic components then estimating the parameters of the single factor should 
lead to much improved results. 


3. It is possible to consider incorporating known information about the future into the 
extrapolation. If, for instance, it becomes known (for data such as the Air Passen- 
ger data) that additional passenger capacity will be provided in Summer months, the 
amplitude of the seasonal demodulates could be increased proportionally. 
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3. Previous Results (1967) 


3.1 Analysis of Buick Sales Data 


The following figures show the results of the application of this procedure to the Buick 
Sales Data time series without trend removal. These are the results originally presented in 
the 1967 paper[11]. It is clear that the trend affected the low frequency demodulates, and 
that the extrapolation is positively biased by the trend. Also, as mentioned above, the plot 
of the demodulates was done on the series before “unwinding” the phase. And, the zero 
frequency demodulate was plotted, then the 7/18 and 7/9 demodulates. This differs from 
the plotted frequecies 7/18, 7/9, and 7/6 plots shown in Figure 2.1.9 (p. 15). 


Original data: gmdata 


Original Data 


2000 


1961 1962 1963 1964 1965 
time 


Figure 3.1: Original Series 
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Figure 3.2: Complex-demodulates (#;,.) to t = 144 
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Initial 1/2 tapering Remodulated Series and Predicted values __ filter length 6, convolved 4 times 
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Figure 3.3: Remodulated (¢ < 144) and predicted (t > 144) Series 
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Figure 3.4: Original and Remodulated-extrapolated Series 
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Initial 1/2 tapering Difference: Original - Remodulated Series filter length 6, convolved 4 times 
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Figure 3.5: Difference: original and remodulated Series 


From inspection of Figures 3.3, 3.4 and 3.5 the predicted values seem to reproduce 
the periodic structure of the series with some indication that the prediction has properly 
adapted to the changing structure of the series. 


DRAFT: 4.1 15 September 2022. 


Conclusion 33 


4. Conclusion 


Based on the analysis using the Air Passenger data, it appears that this method is quite 
effective for extrapolation of series which contain substantial nearly periodic, but time- 
varying, components. The most convincing result is shown in Figure 2.1.7 (p.14). In this 
case the extrapolation began at the start of 1956. The extrapolation for the following 4 years 
was very close to the observed data. Figure 2.1.10 (p. 16) shows that this was partly due to 
the recent stability of the seasonal components. However, when starting the extrapolation 
in 1958 it appeared that the seasonal peak was beginning to be reduced, but that is not 
what happened. 


Comparing these results with the SAS and Spider Financial (p. 17) results which make 
use of the Box-Jenkins ARIMA model shows that at least equal accuracy can be obtained 
without any need for estimation of data-specific parameter values. 


However, highly irregular series such as the Buick Sales Data, cannot accurately be 
predicted. This series was affected by a major strike and other “unexpected” events. Such 
events must detract from the accuracy of any prediction method. 
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5. Appendix 


5.1 Background on the “International Air passenger Data” 


The “International Air passenger Data” time series is frequently used as an example in texts 
on time series and extrapolation. The data are presented in Brown’s book[4] in Appendix 
C, Table C.10 on page 429. The source is stated as: “F.A.A. Handbook of Civil Aviation 
(several annual issues).” A more recent use of the data is in Makidakis/13], pg. 109, where 
other users of the data are referenced. The reported data are a monthly series from January 
1949 through December 1960. The reported units are thousands of passengers. 


The series from which these data were obtained were originally published by the Civil 
Aeronautics Administration (CAA) as the annual Statistical Handbook of Civil Aviation|5]. 
In 1957 the name was changed to the CAA Statistical Handbook of Civil Aviation. Sometime 
between 1959 and 1961 the name was changed to the FAA Statistical Handbook of Civil 
Aviation|6] and the publisher became the Federal Aviation Agency, which, after an number 
of years of conflict, became the successor to the CAA. The last issue to publish these data 
was 1961. In 1962 the monthly data for international travel were discontinued. 


The issues for 1951 and 1952 are missing from the Stanford Library, but the data 
reported in these 2 issues are reported in subsequent issues. 


The first issue to publish monthly passenger data was 1949. The data in this issue run 
from January 1946 through June 1949. After this issue, each issue reported the previous 5 
years of data. 


The data from 1946 through 1960 can be found in the issues of: 


Issue Date Data 
1950: 1946 — 1950, 
1954: 1949 — 1953, 
1957: 1952 — 1956, 
1961: 1956 — 1960. 


Through the 1954 issue, which reported data for 1949 through 1953, the data were reported 
in passenger units. From 1955 the unit was thousands of passengers (three digits). 


It is not made clear in the Handbooks that revisions were made to the data. The 
most common revision was for one previous year. Thus, the data for 1953 as reported in 
1954 were corrected in 1955. It is probable that Brown used the 1954 issue for the data 
for 1946 through 1953. This explains the error for November 1952 and the errors in 1953 
for all months except July. He then used the 1957 issue for 1954 through 1956. This 
explains the errors for the months of January, February, March, August, September, and 
October for 1956. These values were corrected in the 1958 issue. Since 1961 was the last 
monthly data report, there was no chance to correct the 1960 values. Finally, there is an 
unexplained change apparently due to rounding. Rounding to 3 digits started in 1955. The 
value reported in 1954 for December 1950 was 139,532. In 1955 this became 139, but Brown 
reported 140 since he correctly rounded the value from the 1954 issue. 


The largest error was 3 units, or about 1.5%. Frequently, a change of even 0.5% in an 
important economic time series is cause for wide publicity and is viewed as an important 
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indicator of the health of the economy. And, forecasters change their forecasts, Federal 
Reserve action may be taken, etc. The fact that at the next time period the previous 
value is revised by an amount comparable to the original reported change is much less often 
reported, and economists rarely change their minds due to a statistical revision. 


Since this series is now only used as an “example” these errors are of no longer of any 
significance. and, in fact, the extrapolation using the complex-demodulation method is not 
significantly changed by the errors. Still, the fact that in this very popular series, at least 
18 of the 144 values given by Brown were subsequently revised should serve as a reminder 
of the limited accuracy of economic observations. 


Oskar Morgenstern studied the issue of errors in economic data in On the Accuracy of 
Economic Observations [14]). 


5.2 Detailed annotations of the data 


These details are given here since the original publications containing them are not 
widely available. It is curious that the various Public agencies responsible for obtaining 
and reporting economic data have not had sufficient resources to make much use of internet 
technology. 


The Title for the data was always: “Monthly Data — Revenue Passengers Carried 
(Unduplicated).” Early reports gave annual data for both “total” and “revenue” passengers. 
In the 1948 issue there is a footnote reporting that starting in 1947 the carriers were required 
to report only revenue passengers. During the period that both were reported, non-revenue 
passengers were between about 7 and 2 percent of the total. 


The footnotes in each available issue were: 
1949: Source: Bureau of Economic Regulation, CAB. 
1950: Source: Bureau of Air Operations, CAB. 
1953: 4Monthly figures do not add to annual total because of differences in the carriers’ reporting. 
Source: Bureau of Air Operations, CAB. 
1954: ‘Monthly figures do not add to annual total because of differences in the carriers’ reporting. 
Source: Bureau of Air Operations, CAB. 


1955: +4Monthly figures do not add to annual total because of differences in the carriers’ reporting 
and dropping of last 3 digits. 
Source: Office of Carrier Accounts and Statistics, CAB. 

1956: ‘Monthly figures do not add to annual total because of differences in the carriers’ reporting 


and dropping of last 3 digits. 
Source: Office of Carrier Accounts and Statistics, CAB. 
1957: ‘Monthly figures do not add to annual total because of differences in the carriers’ reporting 
and dropping of last 3 digits. 
Source: Office of Carrier Accounts and Statistics, CAB. 


1958: 4Enplaned passengers for 1957 are not strictly comparable to previous years due to changes 
in CAB reporting procedure. 
?Monthly figures do not add to annual total because of differences in the carriers’ reporting 
and dropping of last 3 digits. 
Source: Office of Carrier Accounts and Statistics, CAB. 
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1961: ‘Enplaned passengers for 1957 and subsequent years are not strictly comparable to previous 
years due to changes in CAB reporting procedure. 
?Monthly figures do not add to annual total because of differences in the carriers’ reporting 
and dropping of last 3 digits. 
Source: Office of Carrier Accounts and Statistics, CAB. 


5.3 Example Listing of Airdata in Octave Format: 


comment = "Air Passenger Data. CAA Statistical Handbook" 

numseries = 1; % not used for nonstat_pred 

obsperyr = 12; 

ndemods = [0:24]; 

plotdm = [4 8 12]; % plotted frequencies: 0, pi/4, pi/8 
nfreq = max(size(ndemods)); % number of frequencies for demods 
a = (1.57, -.7]; 

b = [0., .07, .06]; 

nl = 11; % smoothing terms 

nconv = 1; % convolutions 

datlen = 180; 

startex = 84; % extrapolation starts at startex+1 
extlen = 156; % length of extrapolation (from startex+1) 
logt =) a2 

rem_trend = 1; 

di-= La. 

65. 73. 89. 84. 82. 82. 89. 99. 99. 92. 89. 101.... %AIRPAS46 
105. 108. 116. 107. 112. 121. 129. 140. 125. 102. 91. 104.... ‘%%AIRPAS47 
103. 102. 113. 103. 105. 122. 130. 135. 129. 115. 98. 119.... ‘%AIRPAS48 
112. 118. 132. 129. 121. 135. 148. 148. 136. 119. 104. 118.... ‘%AIRPAS49 
115. 126. 141. 135. 125. 149. 170. 170. 158. 133. 114. 140.... ‘%%AIRPAS5O 
145. 150. 178. 163. 172. 178. 199. 199. 184. 162. 146. 166.... %AIRPAS5S1 
171. 180. 193. 181. 183. 218. 230. 242. 209. 191. 173. 194.... ‘AIRPAS52 
195. 194. 234. 232. 230. 242. 264. 271. 236. 209. 178. 200.... ‘%AIRPAS53 
204. 188. 235. 227. 234. 264. 302. 293. 259. 229. 203. 229.... %AIRPAS54 
242. 233. 267. 269. 270. 315. 364. 347. 312. 274. 237. 278.... %AIRPAS5SS 
283. 276. 319. 313. 318. 374. 413. 407. 356. 307. 271. 306.... %AIRPAS56 
315. 301. 356. 348. 355. 422. 465. 467. 404. 347. 305. 336.... %AIRPAS57 
340. 318. 362. 348. 363. 435. 491. 505. 404. 359. 310. 337.... %AIRPAS58 
360. 342. 406. 396. 420. 472. 548. 559. 463. 407. 362. 405.... %AIRPAS59 


417. 391. 419. 461. 472. 535. 622. 606. 508. 461. 390. 432.]; ZAIRPAS6O 
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